function [R,F,A,L,M] = respmtx(ne, dc, sa, ns, sw, x)
% function [R F A L M] = respmtx(ne, dc, sa, ns, sw, k)
%    This function enerates response matrices for slabs in sequence
%    from 1 to ne.  It computes outgoing current responses R, fission
%    rate responses F, absorption rate response A, and leakage rate
%    response L.  Note, is just the outgoing current of vacuum-bounded
%    elements.  Finally, the connectivity matrix M is found via
%    the function connect.

global bcL bcR

R = zeros(2*ne, 2*ne);
F = zeros(1,2*ne);
A = zeros(1,2*ne);

j = 0;

% compute: [ R(11)i  R(21)i;  R(12)i R(22)i ]
% compute: F(left)i F(right)i, A(left)i A(right)i, i=1..ne
k = x(end);
for i = 1:ne

    % B = material buckling
    % B^2 = 0 --> use analytic expression for limit
    %     > 0 --> use trig expressions
    %     < 0 --> use trigh expressions
    
    B2 = (ns(i)/k-sa(i))/dc(i);  % buckling squared
    
    if ( abs(B2) <= 1e-8 ) % buckling nearly zero
        R( (j*2+1), (j*2+1) ) =  sw(i) / ( sw(i) + 4*dc(i) );
        R( (j*2+1), (j*2+2) ) =  4*dc(i) / ( sw(i) + 4*dc(i) );
        R( (j*2+2), (j*2+1) ) =  R( (j*2+1), (j*2+2) );
        R( (j*2+2), (j*2+2) ) =  R( (j*2+1), (j*2+1) );
        F(1,j*2+1) = 2*sw(i)*ns(i);
        F(1,j*2+2) = F(1,j*2+1);
        A(1,j*2+1) = 2*sw(i)*sa(i);
        A(1,j*2+2) = A(1,j*2+1);
        
    elseif ( B2 > 0 )      % positive buckling
        B = sqrt(B2);
        den = 4*dc(i)*B*cos(B*sw(i))+(1-4*dc(i)^2*B2)*sin(B*sw(i));
        R( (j*2+1), (j*2+1) ) = ( sin(B*sw(i))*(1+4*dc(i)^2*B2) )/den;
        R( (j*2+1), (j*2+2) ) = 4*dc(i)*B/den;
        R( (j*2+2), (j*2+1) ) =  R( (j*2+1), (j*2+2) );
        R( (j*2+2), (j*2+2) ) =  R( (j*2+1), (j*2+1) );
        
        den = B*( cos(B*sw(i))*(1+4.*dc(i)^2*B2)+1.-4*dc(i)^2*B2 );
        num = 2*dc(i)*B*(1-cos(B*sw(i)))+sin(B*sw(i));
        F(1,j*2+1) = 4*ns(i)*num/den;
        F(1,j*2+2) = F(1,j*2+1);
        A(1,j*2+1) = 4*sa(i)*num/den;
        A(1,j*2+2) = A(1,j*2+1);
        
    else                   % negative buckling
        B2 = -B2;
        B = sqrt(B2);
        den = 1.0 + 4*dc(i)*B*coth(B*sw(i))+4*dc(i)^2*B2;
        R( (j*2+1), (j*2+1) ) =  (1.0-4*dc(i)^2*B2)/den;
        R( (j*2+1), (j*2+2) ) =  4*dc(i)*B/sinh(B*sw(i))/den;
        R( (j*2+2), (j*2+1) ) =  R( (j*2+1), (j*2+2) );
        R( (j*2+2), (j*2+2) ) =  R( (j*2+1), (j*2+1) );
        % changed -1.-4*dc in den from 1.-4...
        den = B*( cosh(B*sw(i))*(4.*dc(i)^2*B2-1)-1.-4*dc(i)^2*B2 );
        num = 2*dc(i)*B*(cosh(B*sw(i))-1)-sinh(B*sw(i));
        F(1,j*2+1) = 4*ns(i)*num/den;
        F(1,j*2+2) = F(1,j*2+1);
        A(1,j*2+1) = 4*sa(i)*num/den;
        A(1,j*2+2) = A(1,j*2+1);
        
    end
    j = j + 1;
end

L = [0 0];
% leakage is just outgoing currents at boundaries
if bcL == 0
    L(1) = x(1);
end
if bcR == 0
    L(2) = x(end-1);
end

M = connect(ne,bcL,bcR);
R=sparse(R); M=sparse(M);
